Bounds for the Superfluid Fraction from Exact Quantum Monte Carlo Local Densities 
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For solid *He and solid P-H2, using the flow-energy-minimizing one-body phase function and exact 
r = K Monte Carlo calculations of the local density p(r), we have calculated the phase function, 
the velocity proflle and upper bounds for the superfluid fraction fg. At the melting pressure for 
solid ^He we find that fs < 0.20 — 0.21, about ten times what is observed. This strongly indicates 
that the theory for the calculation of these upper bounds needs substantial improvements. 

PACS numbers: 67.80.-s, 67.90.-fz, 67.57.De 



The recent observation of non-classical rotational in- 
ertia (NCRI) effects in bulk solid *He,^ and the experi- 
mental investigation of condensates in optical lattices^ 
have increased theoretical interest in the study of su- 
perfluid properties in presence of translational broken 
symmetry. The interpretation of NCRI effects is still 
controversial and arguments against bulk supersolidity 
in solid ^He have appeared in the literature. To date 
there are microscopic indications of supersolidity from 
exact simulations only in presence of a finite concentra- 
tion of vacanciesi^ The presence of ground state vacancies 
in bulk solid ^He also is controversial)^ Superfluidity of 
interfaces have been also invoked"*, but the quantitatively 
similar presence of NCRI effects in solid ''He confined in 
different porous mediaii^ seems to be difficult to reconcile 
with this single mechanism. 

Here we return to the T = K analysis initiated by 
Leggett in which an upper bound for the superfluid frac- 
tion fs was deduced. In 1970, Leggett studied superflow 
in a solid confined to an annulus3 He considered that, 
under rotation, the ground state wave function ^'o be- 
came Vl* = exp(i</)), with a phase function (jj that was 
a sum of terms depending only upon a single longitu- 
dinal coordinate (what we call a one dimensional one- 
body phase function). He noted that this would yield 
only an upper bound to the T — K superfluid fraction 
fs, but he did not evaluate it. This bound depends on 
the averaged density, p{u) = J d^p(f), where p{r) is the 
local density, u is the longitudinal coordinate and ^ is 
a suitable set of transversal coordinates. By modelling 
p(r) as a sum of Gaussians centered on the lattice sites - 
the so-called Gaussian Model (CM) - Ref. [l3 evaluated 
Leggett's bound. Ref. [lli extended Leggett's bound by 
using the CM and an improved variational ansatz: the 
phase is a function of r, not merely of the longitudinal 
coordinate. The most recent calculation a^^i^^ gave an up- 
per bound for fs on the order of 2%, in good agreement 
with experiments on solid ^He. This bound, however, 
depends very strongly on the value of the width of the 
Gaussians so that this result is model dependent. 

The present work reports calculations of upper bounds 
for fs that, for the first time, use consistent T = K 
densities computed with modern exact Quantum Monte 



Carlo (QMC) techniques. The new results are very dif- 
ferent: the exact QMC one-body density gives an upper 
bound on the order of 20%, for ^He near the melting 
density. This is some ten times higher than observed ex- 
perimentally. Independent of the interpretation of NCRI 
effects in solid ^He, this indicates that the theory for 
calculatiing these upper bounds should be re-examined, 
at least for strongly interacting systems. This judgment 
could be more negative in case the NCRI effects mea- 
sured in solid ''He are not an equilibrium property but 
are induced only by the presence of disorder in the sys- 
tem. 

Ref. 11 showed that a one-body phase function (j){f) 
minimizes the flow energy E = m/2 J drp{r)v1{r), where 
Vs{r) — h/mV(t)(f), when the continuity equation 



W ■ [p{r)vs(r}] = 



(1) 



is satisfied. Thus, on imposing a uniform velocity vq on 
a known p(r), one can obtain Vs(f) and then an upper 
bound for the superfluid density from 



PsV^V < / drp{r)vl{r) 



(2) 



We have computed the local density p{r) in solid ''He 
and in solid P-H2 using an extension of the Path Integral 
Ground State method.^* Speciflcally, we compute p{r) = 

(*o I T.f^i5{r~ n) I «'o), using the exact T = K 
projector Shadow Path Integral Ground State (SPIGS) 
methodi^ In a projector method the exact ground state 
is expressed as the imaginary time evolution of a trial 
variational state; in the SPIGS method this trial state is 
chosen to be a Shadow Wave Function (SWF).— It is 
well documented^ that a SWF gives presently the best 
variational representation of ^He both in the liquid and 
in the solid phase. With R = {fi, .., r/v} representing the 
many-body coordinates, we have^^: 

*o(i?)- lim I dR' G{R,R',t)-^^^^{R!) , (3) 

where G{R,R',t) =< i?| exp(~Ti7)|i?' > is the exact 
imaginary time propagator and 'i>^^^ {R') is the SWF. 



2 



In general, the exact G{R, R' ,t) is not known; how- 
ever, via the path integral representation one can express 
G{R,R',t) as a convolution of many (for example, M) 
short imaginary time propagators G{R, R' , St), where 
St = T /M , and accurate approximations are available for 
G{R, R' , 6t). The imaginary time evolution is stopped 
when r is sufficently large that averages are stable with 
respect to the value of r. This means that '^{R,t) has 
essentially zero overlap with the excited states. In this 
way averages computed with r) for the quantum 

system of N interacting atoms are equivalent to canonical 
averages for a classical system of special interacting open 
polymers of 2M + 1 particles per actual atom, where M 
is the number of projection steps (convolutions) in imagi- 
nary time.^^ Ground state averages of diagonal operators 
can be computed via an efficient Monte Carlo sampling of 
the many-body coordinates for the imaginary time path; 
this is obtained by interpreting the short imaginary time 
propagators and the trial state as probability densities in 
an extended space of N{2M + 1) coordinates. 

In the SPIGS method the imaginary time evolution 
starts from a SWF, which can describe both the liquid 
and the solid phase with the same functional form: spa- 
tial order in the solid phase is produced by spontaneously 
breaking the translational symmetry of the Hamiltonian. 
(To compute the density a special class of moves is made, 
discussed below). With S = {si,..,sn} representing 
the set of subsidiary "shadow" variables Si, we have: 
^^^^{R) = f dS <j)iR)K{R,S)(l},{S). Here 0(i?) con- 
tains the explicit part of the interparticle correlations, 
which are assumed of the Jastrow form (two-body cor- 
relation function); the kernel K{R,S) is often taken as 
a product of gaussians K{R, S) = ni=i e"'^''"'"'*'! , but 
better representations are knowufi^ that correlates each 
coordinate fi with its associated shadow variable s^; and 
(/)s(i?) is another Jastrow factor, which accounts for cor- 
relations between the shadow variables. Integration over 
the subsidiary (shadow) variables introduces, implicitly, 
effective correlations between the particles, which are not 
limited to pair or triplet terms; in principle, terms of 
all order are generated. When the density of the sys- 
tem exceeds the melting density of the solid, such corre- 
lations become so strong that the solid phase stabilizes 
with an energy below that of the (now metastable) liquid 
phasoi^iii; there is no need to introduce ad hoc equilib- 
rium positions to locate the atoms in the solid phase, as 
usually done within the standard variational theory. 

The SWF technique has been shown to be very pow- 
erful but, as in all variational computations, it is intrin- 
sically limited by the functional form of the trial wave 
function. A way to overcome this limitation is to use the 
SWF as a guiding function in an algorithm that converges 
to the exact wave function. This has been obtained with 
the SPIGS method, which inherits the above-mentioned 
properties by projecting out the ground state wave func- 
tion from the imaginary time evolution of a SWF.^^ This 
is especially relevant in the study of a quantum solid, 
where exchange or more complex zero-point phenomena 



could play a relevant role in the determination of prop- 
erties like the local densityi^. Since our wave function 
is translationally invariant the local density turns out to 
be a constant after a very large number of Monte Carlo 
steps due to the fluctuations of the center of mass of 
the system. In order to mimic the behavior of a macro- 
scopic solid held fixed in the laboratory frame we perform 
the SPIGS computation at fixed center of mass: only 
moves that involve all the particles are proposed to the 
Metropolis acceptance test such that when a particle fi 
is moved by S all the others — 1 particles are moved 
by ~6/{N — 1). This is done for each time slice. 

We have computed the local density p{f) in hep solid 
"'He near the melting density (p = 0.029 A""^), at p = 
0.0353 A^'^ and also in hep solid P-H2 at its equilibrium 
density p — 0.026 A~'^, where recent experimental results 
seem to exclude the presence of NCRI effects^i The SWF 
used for the initial state of the imaginary time evolution 
is a standard SWFl&iil. Our SPIGS method employs 
the pair-product approximation^^ for the imaginary time 
propagators. The imaginary time step used in these sim- 
ulations was St = K~^. With AI — 15 time slices this 
gave an imaginary projection time of M5t — 0.1875 K~^, 
and a total of 33 particles in the open polymers (this in- 
cludes two shadow particles). However, only the 11 inte- 
rior time-slices were used to calculate the ground state lo- 
cal density. We have verified that such 6t is accurate for 
solid *He by checking the convergence of diagonal prop- 
erties computed also with St = and K~^. 
The potential used to model the interaction between "^He 
atoms was a standard Aziz potential^^, whereas the po- 
tential of Silvera and Goldman^i was used for P-H2 . 

We first use our QMC densities to evaluate the less ac- 
curate upper bound for fs formulated in Ref. 0, based on 
a one-body phase function that depends on one spatial 
dimension. The averaged local one-dimensional densities 
p{u) along the three axes u = x,y, z oi the hep simu- 
lation box are shown in FigUl they correspond to the 
local density p^f) integrated along the plane perpendic- 
ular to the chosen axis. p(r) has been computed on 64'^ 
real-space points in a cell of dimension (a, \/3a, -\/8/3a). 
There are 45 equivalent cells of this kind in the hep sim- 
ulation box with N=180 atoms; averaging the local den- 
sities in each cell has been used in order to improve the 
statistics. FigH] shows that, of the three /9(u)'s, p{z) (for 
the direction perpendicular to the basal plane) has the 
strongest oscillations around the average density p in the 
system. This method of averaging includes peaks from 
different planes; for p{y), the peaks from different planes 
overlap. As seen in FiglU near melting {p — 0.029 A"'^) 
the peaks are so broad that the average gives a single 
broadened "peak" in p{y). However, at higher pressure, 
where the atoms are more localized {p = 0.0353 A'^), 
the contributions of the individual planes can be partially 
resolved. For each of the three cases depicted in FiglU 
we have computed the upper bounds for flow along x, y, 
and z. In each case, the lowest of these upper bounds 
is associated with flow along z, because the large oscilla- 
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FIG. 1: Averaged local densities along the axis of the simulation box for hep solid ''He at p = 0.029 (solid line), hep solid 
''He at p = 0.0353 (dotted line), hop solid P-H2 (dot-dashed line) at p = 0.026 A~^. The basal plane is perpendicular to 
the z direction; the x axis and y axis are respectively parallel to the FK and FM direction. The choice of the origin along the 
axis is arbitrary. The periodicity of the oscillations in the averaged local densities depends on the lattice constant and changes 
with the average density p. 



tions in p{z) give deeper minima. These lowest Id upper 
bounds are fs < 0.384 for hep solid "He at p = 0.029 
A^'^ (flow along X gives fs < 0.939 and flow along y gives 
fs < 0.799), fs < 0.164 for hep solid "He at p = 0.0353 
A^'' (flow along X gives fs < 0.839 and flow along y 
gives fs < 0.648) and fs < 0.054 for hep solid P-H2 at 
p = 0.026 A~^ (flow along x gives fs < 0.686 and flow 
along y gives fs < 0.458). 

An improved upper bound for fs is obtained by using a 
one-body phase function that depends on all three spatial 
dimensions and the local density p(r) given by SPIGS. 
Note that the one-body phase is the sum of an imposed 
uniform flow and a one-body backflow that varies on the 
atomic scale as follows from Eq.(IT]). Eqs.(IT][2|) are solved 
in Fourier space. We have studied fs for the hep lattice, 
as actually realized experimentally. In this case the crys- 
tal structure is a Bravais lattice with a basis of two. By 
placing the origin midway between the basis atoms, one 
may employ only real fourier transforms. The results are 
given m Table [11 for hep "He at p = 0.029 A^^^ for hep 
"He at p = 0.0353 A"^ and for hep p-Hs at p = 0.026 
A""^; the number of fourier components associated with 
each reciprocal lattice vector basis is given by 2P + 1 . In 
the table the accuracy of the QMC density does not war- 
rant four decimal places, but the same QMC density was 
used for each calculation, so that relative values have four 
place accuracy. The fs shown in Table |T] have been ob- 
tained from p{r) computed with an imaginary projection 
time of M6t = 0.1875 K"-'^; we have checked the conver- 
gence of these results by computing the local density with 
different imaginary projection times: AISt — 0.15 K^^ 
and M5t = 0.2625 K~'. In Table U we have shown also 
fs computed with SWF; the variational fs turn out to be 
always lower as a consequence of the larger degree of local 
order—. This is consistent with the larger Lindemann's 
ratio 7 = 0.257(4) obtained with SPIGS at p = 0.029 
A~'^ as opposed to the one, 7 = 0.242(2), obtained with 
SWF; the experimental value is 7*=^^ = 0.263(6)^^2 Note 



that, to within numerical accuracy, all three systems are 
isotropic. Such isotropy of fs need not hold for a hexag- 
onal lattice, but near-isotropy is consistent with a Gaus- 
sian Model (GM) for the hexagonal latticeJ^ Although 
inspection of the tables indicates that these three systems 
are nearly isotropic, in each case visual inspection of the 
data (connected by straight Hues) indicates that /s(OOl), 
associated with flow along z, is extrapolating to a value 
a few percent lower than for the other two. Similarly, vi- 
sual inspection indicates that /^(lOO) and /^(OlO) appear 
to extrapolate to the same value, indicating isotropicity 
for flow along the basal hep plane. 

Converged values for fs for "He at melting density p = 
0.029 A-3 (/,_,, = 0.220, fs^yy = 0.218, /,,^, = 0.211) lie 
in the range fs = 0.21 - 0.22. In terms of the GMi^, 
this corresponds to an rms width of ti ~ 0.1485 d, where 
d is the nearest-neighbor distance for the hep lattice. 
Note that in the GM, a = 0.1414 d gives fs = 0.133 and 
cr ~ 0.159 d gives fs = 0.303, so that fs is a very sensitive 
function of the localization. Variations in fs according to 
direction are likely due to numerical uncertainties in the 
QMC density p(r), of the order of 3% at this density. 
This uncertainty is unlikely to be the source of the factor 
of ten difference between the one-body calculation de- 
scribed above, and the experimental values on the order 
of 1-2%,' The system is therefore less localized than pre- 
vious calculations had indicated)^ and this is reflected in 
the larger value of fs , which increases with the extent of 
delocalization, in principle to a maximum of unity (the 
fluid state). A possible source of discrepancy lies also 
in the usage of the GM in the previous calculationSf^^ 
this model in fact loses accuracy when we consider the 
regions of the minima of p(f)', here deviations from the 
exact QMC local densities greater than 100% are found. 
However these discrepancies seem to have only a small 
effect on the derived upper bounds fs'. by fitting the 
QMC averaged local densities with the GM we obtain 
for "He at p = 0.029 A'^ a rms width a = 0.1486 d 
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TABLE I: Superfluid fraction upper bound, computed with SPIGS and SWF, for hep solid *He and hep solid P-H2 at different 
densities for varying basis sizes (2P + 1)^ in fourier space. 
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1 


0.4055 


0.3801 


0.5694 


0.5015 


0.2517 


0.2253 


0.3968 


0.3736 


0.1553 


0.1336 


0.2653 


0.2238 


2 


0.2884 


0.2733 


0.3340 


0.2590 


0.1267 


0.1127 


0.1643 


0.1457 


0.0526 


0.0436 


0.0736 


0.0512 


3 


0.2464 


0.2389 


0.2618 


0.1852 


0.0860 


0.0792 


0.1004 


0.0846 


0.0263 


0.0225 


0.0332 


0.0195 


4 


0.2310 


0.2267 


0.2346 


0.1573 


0.0705 


0.0667 


0.0770 


0.0626 


0.0175 


0.0155 


0.0207 


0.0106 


5 


0.2249 


0.2219 


0.2224 


0.1444 


0.0637 


0.0614 


0.0662 


0.0526 


0.0137 


0.0126 


0.0154 


0.0071 


6 


0.2224 


0.2199 


0.2166 


0.1379 


0.0605 


0.0589 


0.0606 


0.0473 


0.0119 


0.0112 


0.0127 


0.0054 


7 


0.2213 


0.2190 


0.2137 


0.1345 


0.0588 


0.0576 


0.0576 


0.0444 


0.0109 


0.0105 


0.0112 


0.0045 


8 


0.2208 


0.2186 


0.2124 


0.1327 


0.0580 


0.0569 


0.0558 


0.0427 


0.0104 


0.0101 


0.0103 


0.0039 



and at p = 0.0353 A'^ a rms width a 0.1274 d; in 
the GM^^'^'^ these rms widths correspond respectively to 
fs ~ 0.22 and fs « 0.05, in good agreement with the 
results in Tabled SoHd ''He at p = 0.0353 A'^ and soHd 
P-H2 at p = 0.026 A-3 are more localized systems, and 
give correspondingly lower values for fs (see Table 
however even if there seems to be no indication of NCRI 
effects in solid p-H2^ the upper bound is still over 1%. 
We have also computed the upper bound for the super- 
fluid fraction for ^He on an fee lattice; at 0.029 A~'^ the 
results are close to those obtained for the hep lattice. 

To conclude, we have performed calculations of the up- 
per bound for the superfluid fraction using an optimized 
one-body phase function and state-of-the-art QMC one- 
body densities. We find an upper bound for the super- 
fluid fraction of 20% for solid ''He at density of 0.029 
A~'^ (near the melting pressure), significantly higher than 
what has been observed. Therefore the usage of exact lo- 
cal densities has definitely shown that the accuracy of the 



upper bounds for the superfluid fraction is of little util- 
ity in the present form, at least for a strongly interacting 
quantum system like solid ^He and solid P-H2; in the case 
of condensates in optical lattice the conclusions could be 
different. However, given an additional degree of freedom 
by permitting a many-body phase function, it should be 
possible to lower the flow energy for the case of the solid, 
and therefore to lower the corresponding superfluid frac- 
tion, perhaps by a considerable amount. We have devel- 
oped such a theory^^ for a two-body phase function, and 
the results should also be applicable to condensates in 
optical lattices and to disordered solids. 
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